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Abstract. In this paper we introduce a modified lattice Boltzmann model (LBM) with the capability 
of mimicking a fluid system with dynamic heterogeneities. The physical system is modeled as a one- 
dimensional fluid, interacting with finite-lifetime moving obstacles. Fluid motion is described by a lattice 
Boltzmann equation and obstacles are randomly distributed semi-permeable barriers which constrain the 
motion of the fluid particles. After a lifetime delay, obstacles move to new random positions. It is found 
that the non-linearly coupled dynamics of the fluid and obstacles produces heterogeneous patterns in fluid 
density and non-exponential relaxation of two-time autocorrelation function. 

PACS. 47.11. +j Computational methods in fluid dynamics - 05.70.Ln Nonequilibrium and irreversible 
thermodynamics 

1 Introduction the temperature is decreased (or density increased) 

Within this picture, long-time relaxation is often associ- 
Slow relaxation to local equilibrium is a hallmark of com- ated with thfi appearance of sp ace-time heterogeneities, 
plex system behaviour, with many examples in physics, which can be tentatively attributed to self-trapping ef- 
material science, and biology I, From a many-body point i {J molecules get trap ped into 'cages' formed by 

of view, the emergence of long-time relaxation seems to be high . density aggregations of other groups molecules. While 
related to the gradual confinement of the system in lower- the landscape picture necesS arily calls for many-body in- 
dimensional regions ('slow' manifolds) of phase-space, sur- ves tig at ions (molecular dynamics, Monte Carlo and 
rounded by isolating barriers of increasing amplitude as 
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various types of lattice-'glasses' [3], for the case of glasses), reads as follows (time-step made unit for simplicity): 
dynamic heterogeneities can also be interpreted in terms 

Mr, t) - Mr - Ci , t - 1) = -w [fi - ff ] (r - c h t - 1) (1) 

of mutual interactions between fluid molecules of differ- 

iPi T£ u i • i • ii where fi{r,t) = /(r,v = Cj,t) is a discrete distribution 

ent mobilities 5^ If such a picture is accepted, it is then 

... function of particles moving along the direction i with 

natural to attempt a description m terms ol much sim- 

, rr . , , , , n i i\ i mi • • discrete speed Cj. The right-hand side represents the re- 

pler effective single- body (mean-held) approaches. I his is 

. laxation to a local equilibrium /? in a time lapse of the 

precisely the conceptual framework of the present work. 

We develop a mesoscopic model of dynamically het- or< ^ er w 

_ . , , . , . , , . _ . . The equilibrium distribution functions ff are expressed 

erogeneous fluids based on the lattice Boltzmann equation 

n nx^\ x • • • i i • , • i -i • i as series expansions of the Maxwellian distribution func- 

(LBE). LBE is a minimal kinetic equation describing styl- 

. . . . , . . i , • i tion, up to second order with respect to the local velocity 

lzcd pseudo-molecules evolving m a regular lattice accord- 
ing to a simple local dynamics including free-streaming, ^ 

collisional relaxation and (effective) intcrmolccular inter- f e = pwAl H H v - J J—L\ (o) 

actions jSJ- LBE has proven extremely successful for a 

The local density p — p(r,t), as well as the local velocity 

variety of complex flows, but its applicability to disor- 

u which enter eq. (J2J) , are calculated from the distribution 

dered fluids with long-time relaxation appears more prob- 

functions as follows: 

lematic 0. A crucial ingredient to produce dynamic het- 
erogeneities is (dynamic) geometrical frustration, i.e. the ^ ~ ^ 

i 

introduction of constraints which reduce the phase-space u _ _ j\ c . 

^ i 

available to the fluid system. To model these effects, the 

Wi is a set of weights normalized to unity, c s the lattice 

dynamics of the LBE molecules includes an interaction 

sound speed, and finally, / stands for the unit tensor. In 

with space-time dependent obstacles hindering their mo- 

the following, we shall refer to a one-dimensional lattice 
bility. As we shall see, the coexistence of 'slow' (the ob- _ 

of size L with = —1,0 + 1 and c s = 

stacles) and 'fast' (LBE molecules) makes the dynamics 



of the present model depart significantly from simple fluid 
behavior. 



We generalize the above eq. Q as follows 



fi(r,t)- [l-p(r +f>*-l)] fi(r,t-l) 
p(r-^,t-l)\Mr-c h t-l) 



2 

2 The model -u [h - ft) (r - a, t - 1)] (5) 

In the present paper, we use a modified lattice Boltzmann The variables p(r ± t) live on the lattice links and vary 
equation. Standard LBE with a single relaxation time (Hj in space and time (non-quenched disorder |1()|L We choose 
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Fig. 1. Schematic representation of dynamical role of obstacles 
during propagation of distribution functions. 

them in the form of binary fields taking only the values 
p = 1 and p = pt with < pt < 1. Links with p = p t 
act as semi-permeable obstacles which control the prop- 
agation rate between neighboring sites in a way which 
preserves the total density conservation. Figure Q illus- 
trates in a pictorial way the role of obstacles. Obstacles 
are characterized by their permeability, p t , and concentra- 
tion, c = O/L, O being the number of obstacles, which is 
kept fixed in time. 

The dynamics of the obstacles is the following. At time 
t = 0, the number of obstacles, O, is chosen and their 
positions are randomly selected. At subsequent times t n , 
n = 1,2,.., where 

tl = T 

tn+1 = tn + T(t n ) (6) 

r(t„) = m^Toe^fed™^)- 2 ! -1 )^ 

the obstacles are shifted rightward to new randomly cho- 
sen positions r i (i„ +1 ) = r,(i„) + r' j (t n+i ), j = 1,...,0, 
where r'j(t n -\-i) = int[s) + 1, s being a random number 
drawn from an uniform distribution in the range [|, 4^] 
and d — jy is the mean inter-obstacle distance. In the 
above, int(s) denotes the integer part of the variable s. 
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The motion of each obstacle can be seen as a continuous 
time random walk 11]. In eqs. JBJl To is the first lifetime 
and the quantity 

m (tn) = {PmaxiUi) ~ Pmin{t n )) / p Q (7) 

is the relevant order parameter. Here, Pmax(t n ) and p m in(t n ) 
are the maximum and minimum values of fluid density at 
time t n , respectively, and po is the average fluid density. 
Non zero values of m are the prime indicators of depar- 
ture from ideal fluid behavior (see below). The quantity 
To has to be generally much larger than the time scale of 
fluid motion so as to characterize obstacles as the slowly 
moving species of particles. The dependence of t on the 
order parameter m is intended to slow down the dynamics 
of obstacles in the presence of density contrasts (non zero 
values of m). Indeed, in the expression of r(t n ) in © it ap- 
pears a singularity at m = 2. The reason is the following. 
We assumed that, on the average, p m ax — Pa — Po — pmin 
and later numerically verified this assumption to be cor- 
rect. Since it must be p m in > 0, this means that in the 
limit p m i n — > (pmax — * 2/?o), namely m — > 2, the lifetime 
of the obstacles must diverge, r(t n ) — > +oo. 

In moving obstacles we have taken into account the pe- 
riodicity of the lattice. When a link is already occupied, 
the moving obstacle is shifted to the nearest neighbor link. 
The fact that obstacles move rightwards has no effect on 
the overall fluid dynamics, which depends only on the po- 
sition of obstacles, and not on their motion. The equation 
of motion of the obstacles is: 

p{r + r',t n+1 )=p{r,t n ), n=l,2,... (8) 
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where the distance r' has been previously defined. The 
eq. © is coupled to eq. JSJ) via r(i n ), which depends on 
density through the order parameter m = m(p). 

Before proceeding further, a few distinguished limits 
in the phase-plane p t — c are worth a brief comment: 

- Fluid limit: Along the lines c = (obstacle-free) and 
p t — 1 (transparent obstacles) the system behaves like a 
standard LBE fluid. 

- Slow- fluid limit: Along the line c — 1, the system be- 
haves like a standard LB fluid, only with a rescaled speed 
Ptu. Thus, as pt — > 0, this slow fluid goes smoothly into 
the 'frozen' limit represented by the corner (p t = 0, c = 1), 
where the system is frozen to its initial configuration. 

- Arrest limit: Along the line p t — the system is li- 
able to develop density accumulations in a finite time at 
the locations where the obstacles reside. The actual onset 
of these density pile-ups depends on the lifetime of the 
traps, as well as on their concentration. It is in fact clear 
that dilute obstacles with short lifetimes permit the sys- 
tem to release density pile-ups, thereby avoiding strong 
density accumulations. In actual practice, an 'arrest line' 
of the form ca = CA{pt, To) marks the upper value of the 
concentration leading to configurations without accumu- 
lation. Intuitively, ca is an increasing function of pt and 
to- 

The macroscopic limit of the present model remains an 
open question mainly because the obstacle populations is 
generally smooth, so that a standard Chapman-Enskog 
analysis is difficult to apply. Here we only present a few 
remarks which apply to the limit pt — > 1. To leading order 
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in the parameter l—p t , the main features of density pro- 
files can be explained in terms of the following modified 
continuity equation: 

d t p + d a (ppu) = 0. (9) 

By writing the above as: 

d t p + d a {pu) = -pud a p - (p - l)d a {pu) (10) 

we recognize two extra compressibility terms on the right 
hand side. The limit p — 1 annihilates both extra com- 
pressibility contributions, as it must be, since this is the 
standard LBE situation. The case p = p t everywhere (c = 
1), leads again to a standard LBE, only with a rescaled 
speed u — > p t u, hence no effective extra-compressibility 
Genuinely extra-compressibility is therefore associated to 
spatial changes of the permeability field p(r, t) : d a p ^ 
(this derivative must be intended in the sense of distribu- 
tions, since p{x,t) is generally not smooth). 

3 Numerical results 

We performed a series of numerical simulations on a lat- 
tice with L = 1024 grid-points and to — 1.5. We found 
that results are not dependent on u), which was varied in 
the range [1,1.7], corresponding to a local collisional re- 
laxation timescale r c ~ 1/lo approximately in the range 
[0.58, 1] in time step units. The parameter To was fixed at 
To = 50. The initial condition is p{r) = po (l + £(r)) where 
£(r) is a random perturbation uniformly distributed in the 
range [—0.1 : 0.1] and po = 1. 

As a first task, we determine the phase-diagram of our 
model in the pt — c plane. We spanned the pt — c plane and 
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Fig. 2. Left panel: The order parameter m as a function of pt for different concentrations c = 0.5 (•), 0.78 (o). The straight 
line has slope —3/2. Right panel: The order parameter m as a function of c for different values of pt = 0.7(«), 0.9(o). 



found that departures from ideal fluid behavior (m =/= 0) 
are observed for every value of (pt, c) with p m i n < Pt < 1 
and < c < 1. The limiting value (arrest value) p m in indi- 
cates the lowest permeability, below which the simulation 
is disrupted by excessive density pile-up. 

As expected, the arrest value p m i n decreases with de- 
creasing c: We found p m in — 0.52 for c = 0.5 and p m in — 
0.61 for c = 0.78. The behaviors of m as a function of 
Pt and c are shown in fig. [21 In the vicinity of p m im we 
observe a power-law behavior m ~ (p t — p m in)~ a with 
a — 3/2. Away from p mm , m grows like {p max {c) - Pt) b , 
with b ~ 1/2, which is essentially a mean field theory ex- 
ponent, and Pmax{c) — > 1 as c —> 1. The dependence of 
the order parameter on c at fixed p t is non-monotonic, 
with a kink around c ~ 0.02 and a sharp collapse towards 
the 'slow-fluid' limit, c = 1, which appears to be reached 
in a highly non-perturbative way. To be noted that even 
a single obstacle c = 1/L ~ 0.001 is sufficient to gener- 
ate sizeable non-zero values of the order parameter m. In 
summary, these results indicate that fluid behavior is re- 



covered smoothly, but shows very small resilience towards 
non-zero values of the obstacle concentration c and imper- 
meability 1 — pt ■ 

We also considered the dependence of the order pa- 
rameter m on the lifetime To. Figure [3] shows the behavior 
of m as a function of the parameter tq in the case with 
c = 0.78, pt = 0.7. It appears that m decays with in- 
creasing To to reach a constant value m(oo). The reason 
of this behavior is the following. In the limit of static dis- 
order (to — ► oo) density profiles are flat among obstacles 
and show sharp contrasts across them. Density accumu- 
lations are prevented in the long time limit by the flow 
which flattens density among obstacles. When reducing 
To it is more difficult for the flow to smooth density pro- 
files and the motion of obstacles produces larger density 
contrasts in the system as obstacles are shifted thereby 
causing m to increase. In the limit of rapid motion of ob- 
stacles (to — ► 0) it happens that density contrasts are so 
steep that m diverges. Further simulations with p t in the 
range [0.7, 0.9] indicate that that m(oo) is close to the 
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Fig. 3. The order parameter m as a function of to for c = 0.78 
and p t — 0.7. The full line is a guide to the eye and represents 
eq. |TTJ. 

value (1 —p t )/((l +pt)/2), the latter being an estimate of 
the density variation across a single obstacle based on the 
continuity equation J5J). It is interesting to observe that 
the dynamic component of the order parameter m obeys 
a scaling law of the form 

m(r ) ~ ra(oo) = a 

m(oo) ( T0 _ T(W )i/ 2 ^ > 

again with a mean-field like exponent 1/2, a and To crit 
being two parameters which depend on p t . 

We next turn to the analysis of the dynamics of the 
model, notably through a study of the time and space 
correlation functions. To this purpose, we choose c = 0.78, 
p t — 0.7 and tq = 50 in order to analyse a region of 
the phase diagram with clear departures from ideal fluid 
behavior (to ~ 0.75, see figs. I2I3JI . 

We verified that the resulting motion of obstacles can 
be described in terms of a directed random walk. Indeed, 
we have computed the time evolution of the distance D 
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Fig. 4. The distance D travelled by obstacles as a function of 
times t„. The straight line has slope 1. 
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Fig. 5. The density autocorrelation function h(t, t w ) as a func- 
tion of t for different waiting times t w — 1(A), 52(o), 104(*), 
157(o), 211(*). h(t,t w ) is shown also in the case without ob- 
stacles for t w = 211(«). 

travelled by obstacles defined as 

y° 1 rJt n ) -r,-(Q) 
D(t n ) = ^i. (12) 

The quantity D is plotted in fig. 0] as a function of times 
t n . From this figure a power law D ~ t is well visible. 
The exponent 1, characterizing the time behavior of D, is 
consistent with the fact that the motion of obstacles can 
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be seen as a directed random walk for which the exponent over the whole lattice. The corresponding probability dis- 

is known to be 1. This result confirms that the obstacle tribution functions (PDF) are shown in fig. In the ini- 

motion is independent on fluid coupling and on obstacle tial stage the PDF quickly forgets the initial distribution 

shifting. and develops two peaks around the initial minimum and 

We focused our attention on the two-time density au- maximum densities 0.9 and 1.1 respectively. The appar- 

tocorrelation function defined as ent phase separation at time i — 50 is just due to the 



w. , n _ < 5p(r,t w )Sp(r,t w + t) > r . . effect of obstacles which cause density to accumulate on 

( ' W ' ~ <p(r,t w )* > r 



where Sp(r) = p(r) — po is the density fluctuation around 
its spatial average, < ... > r denotes sum over the whole 
system and t > t w . Multiplying two configurations (one 
at time t = t w , the other at t w + t) together site by site 
is sufficient to produce a satisfactory average of h(t,t w ) 
because, due to the large number of sites [L — 1024), the 
above procedure is equivalent to obtain h(t,t w ) from an 
ensemble average of the system \12\ . 

The function h(t,t w ) is not normalized to unity at 
t = due to the definition (eq. (|13[0. This definition is 
common for correlations of fluctuations (dp) since small 
values of fluctuations at the denominator can give di- 



some lattice sites. Indeed, as soon as the obstacles are re- 
leased this effect tends to fade away. As time unfolds, these 
peaks 'diffuse' with the twofold result of filling up the gap 
at p = po and populating both high and low density tails. 
It is interesting to note that the PDF generated in these 
processes are not Gaussian, hinting at the presence of (dy- 
namic) heterogeneities in the system. In the case without 
obstacles density p gets constant at value po (see fig. at 
time t — 5000) giving a single peak at p = po in the PDF 
(see fig. □ at time t = 5000). 

As a further observable, in fig. [S] we plot the density 
correlation function: 



verging contributions to the correlation function itself. q(x) — — — -. (14) 

< p(r) 2 >r 



The plot for different waiting times is presented in fig. |SJ 
From this figure, evidence of non-exponential relaxation 
is clearly seen, hinting at the presence of long-lived dy- 
namical structures. Also, the dependence on the wait- 
ing time t w indicates the existence of an underlying non- 
equilibrium, non-stationary process. In the case without 
obstacles (standard LBE with p = 1 in eq. (J5J) h(t, t w ) is 
approximately zero (see fig.|SJ), as it should be. 

Plots of density are shown in fig.|Hlin the range 1 — 500 
for a better view. The qualitative behavior is the same S p (x,t) —< \dp(r, x, t)\ p > r , p= 1,2, ..,5 (15) 



From these pictures, the onset of long-range order is well 
appreciated. In the case without obstacles g(x) is approx- 
imately zero (see fig- El a t time t — 5000). 

In order to gain further insights into the dynamics of 
the density field, we performed a scale-dependent analysis. 
Let dp(r, x, t) = p(r + x, t) — p(r, t) be the density incre- 
ment at a scale x at time t, and define the corresponding 
structure functions as: 
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Fig. 6. Density profiles at different times on a portion of the lattice. At time t = 5000 the profile is shown also in the case 
without obstacles (thick line). 
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Fig. 7. The probability distribution functions (PDF) of p at different times. At time t = 5000 the PDF is shown also in the 
case without obstacles (—•—). 



0.015 



-.0.01 



_><_ 

CO 



0.005 



-0.005 




0.015 



-0.01 - 



0.005 



-0.005 




0.015 



„0.01 

_><_ 

CD 

0.005 



-0.005 



t=5000 




V- 



XDCOOOOOOOOOOCOGOOOOOO: :: : : — : 



oooooSo***M' 



10 20 30 40 10 20 30 40 10 20 30 40 

xxx 

Fig. 8. The density correlation function g(x) as a function of x for different times. At time t = 5000 g(x) is shown also in the 
case without obstacles (o). 
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Fig. 9. The structure factors S p (x,t) as a function of Ss(x,t) for different values of p = 1(A), 2(o), 3(*), 4(«), 5(*) at two 
times. The straight line have slopes 5/p, p = 1,2, .., 5 from right to left. 



as well as scaling exponents a p via: 



S p {x,t) ~ S p (x ,t)( — ) a ?, 1<X «L (16) 
x 



where xq = int(d). It is well known that a linear depen- 
dence a p = Kp indicates a fractal process of order K 
(K = 1 for smooth, differentiable processes, K = 1/2 for 
standard diffusion), whereas a non-linear dependence onp 
would signal a multi- fractal process instead J3j . Our data 
show that S p (x,i) ~ x Ktyt ^ p , where K(t) is a time-varying 
but scale-independent parameter. This suggests that the 
density diffusion process is a fractal of dimension K(t), 
with K raising in the course of the evolution from 0.08 
to the steady value 0.15. This also hints at a hierarchi- 
cal organization of density peaks and dips, which emerges 
spontaneously from the non-linear coupling between the 
fluid motion and the obstacles dynamics. 

The fractal character is more neatly evidenced by using 
extended self-similarity ^3]- Indeed, in fig. [5] we show the 
plots of S p (x, t) vs S§(x, i) at two times for p = 1,2,.., 5. 
Data points fall on the straight lines of slope 5/p, as ex- 



pected in the case of fractal behavior. The low values 
of K (t) indicate that the density redistribution is a sub- 
diffusive process, i.e. it proceeds more slowly than Brow- 
nian motion |KS| , as it is expected for motion in a hetero- 
geneous, disordered background. 



4 Conclusions 

In summary, we have presented a mesoscopic model of 
fluid motion with random dynamical constraints. Fluid 
motion is based on a lattice Boltzmann model, whereas 
dynamical constraints arc implemented via a control field 
acting as a penetrable barrier for particle motion along the 
links. Despite its simplicity, our model displays some fea- 
tures of disordered fluids, namely: i) Onset of non-zero or- 
der parameter (density contrast) with very small resilience 
to non-zero obstacle concentration c and impermeability 
1— pt] ii) Non-exponential, non-stationary time relaxation 
of density correlation functions; iii) Long-range spatial or- 
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der, with non-Gaussian PDF's of the density distribution, 
iv) Sub-diffusive dynamics of the density distribution. 

Although no strict correspondence with any specific 
physical system can be claimed at this stage, the natural 
physical target of our model are fluids in porous media 
with a dynamic morphology nonlinearly coupled to fluid 
motion. For the future, we plan to investigate the behav- 
ior of the present model for two and three-dimensional flu- 
ids. In the long-term, we would also like to model glassy 
materials, although this is certainly going to require new 
substantial extensions of the simple model presented in 
this work. 

Illuminating discussions with K. Binder, E. Marinari, G. Parisi 
and F. Sciortino are kindly acknowledged. 
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